Snakemake Report

Loading Snakemake Report...

Please enable Javascript in your browser to see this report.

Loading 1.4 MB. For large reports, this can take a while.

Workflow

Click the nodes to obtain details about each step.

Other

Statistics

If the workflow has been executed in cluster/cloud, runtimes include the waiting time in the queue.

Configuration

File Code

method: XR

sample: 
  - "HXA64A1_ATCACG_sub10k"

layout:
  - "single-end"

adaptor_se: "-a TGGAATTCTCGGGTGCCAAGGAACTCCAGTNNNNNNACGATCTCGTATGCCGTCTTCTGCTTG"

adaptor_pe: "-a TGGAATTCTCGGGTGCCAAGGAACTCCAGTNNNNNNACGATCTCGTATGCCGTCTTCTGCTTG -A TGGAATTCTCGGGTGCCAAGGAACTCCAGTNNNNNNACGATCTCGTATGCCGTCTTCTGCTTG"

cutadapt_se: ""

cutadapt_pe: ""

bowtie2: ""

samtools_q_trim_se: "-q 20"

samtools_q_trim_pe: "-q 20 -bf 0x2"

filter: "'^([1-9]|1[0-9]|2[0-2]|X)'"

# Genome parameters 

build: "GRCh37"

species: "homo_sapiens"

datatype: "dna"

release: "102"

#######

sample: 
  - "HXA64A1_ATCACG_sub10k"


build: "GRCh38"

Jobs 1

Input files

resources/samples/{samples}.fastq.gz

Output files

results/{samples}/{samples}.html
results/{samples}/{samples}_fastqc.zip

Conda software stack

fastqc ==0.11.9

Code

"""Snakemake wrapper for fastqc."""

__author__ = "Julian de Ruiter"
__copyright__ = "Copyright 2017, Julian de Ruiter"
__email__ = "julianderuiter@gmail.com"
__license__ = "MIT"


from os import path
from tempfile import TemporaryDirectory

from snakemake.shell import shell

log = snakemake.log_fmt_shell(stdout=False, stderr=True)


def basename_without_ext(file_path):
    """Returns basename of file path, without the file extension."""

    base = path.basename(file_path)

    split_ind = 2 if base.endswith(".fastq.gz") else 1
    base = ".".join(base.split(".")[:-split_ind])

    return base


# Run fastqc, since there can be race conditions if multiple jobs
# use the same fastqc dir, we create a temp dir.
with TemporaryDirectory() as tempdir:
    shell(
        "fastqc {snakemake.params} --quiet -t {snakemake.threads} "
        "--outdir {tempdir:q} {snakemake.input[0]:q}"
        " {log:q}"
    )

    # Move outputs into proper position.
    output_base = basename_without_ext(snakemake.input[0])
    html_path = path.join(tempdir, output_base + "_fastqc.html")
    zip_path = path.join(tempdir, output_base + "_fastqc.zip")

    if snakemake.output.html != html_path:
        shell("mv {html_path:q} {snakemake.output.html:q}")

    if snakemake.output.zip != zip_path:
        shell("mv {zip_path:q} {snakemake.output.zip:q}")

Jobs 1

Input files

results/{samples}/{samples}_{build}_sorted_dinucleotideTable.txt
results/{samples}/{samples}_{build}_sorted_nucleotideTable.txt

Output files

results/{samples}/{samples}_{build}_sorted_dinucleotideTable.png
results/{samples}/{samples}_{build}_sorted_nucleotideTable.png

Conda software stack

conda-forge::r-base ==3.5.1
conda-forge::r-argparser ==0.4
conda-forge::r-ggplot2 ==3.0.0
conda-forge::r-dplyr ==0.8.0.1
conda-forge::r-tidyr ==0.8.2

Code

        Rscript workflow/scripts/nucleotidePlot.r         -i {input.dinuc}         -k 2         -s {params.name}         -f {params.motif}         -o {output.dinuc} 2>{log}

        Rscript workflow/scripts/nucleotidePlot.r         -i {input.nuc}         -k 1         -s {params.name}         -o {output.nuc} 2>>{log}

Jobs 1

Input files

results/{samples}/{samples}_{build}_lengthMode.fa

Output files

results/{samples}/{samples}_{build}_sorted_dinucleotideTable.txt
results/{samples}/{samples}_{build}_sorted_nucleotideTable.txt

Code

        workflow/scripts/fa2kmerAbundanceTable.py         -i {input}         -k 2         -o {output.dinuc} 2>{log}

        workflow/scripts/fa2kmerAbundanceTable.py         -i {input}         -k 1         -o {output.nuc} 2>>{log}

Jobs 1

Input files

results/{samples}/{samples}_{build}_lengthMode.bed
resources/ref_genomes/{build}/genome_{build}.fa

Output files

results/{samples}/{samples}_{build}_lengthMode.fa

Conda software stack

bedtools =2.29.0

Code

1 2	`bedtools getfasta -fi {input.genome} -bed {input.bed} -fo {output} -s`

Jobs 1

Input files

results/{samples}/{samples}_{build}_sorted_chr.bed
results/{samples}/{samples}_{build}_length_distribution.txt

Output files

results/{samples}/{samples}_{build}_lengthMode.bed

Code

        length="$(awk -v m=0 '{{if(m<$2){{m=$2;l=$1}}}}END{{print l}}'         {input.ld})"

        awk -v num="$length" '{{ if ($3-$2 == num) {{ print }} }}' {input.bed}         > {output}

Jobs 1

Input files

results/{samples}/{samples}_{build}_.bed

Output files

results/{samples}/{samples}_{build}_sorted_chr.bed

Code

1
2
3

        sort -u -k1,1 -k2,2n -k3,3n {input} |
        egrep {params.filt} > {output}

Jobs 1

Input files

results/{samples}/{samples}_cutadapt_se_{build}.bam

Output files

results/{samples}/{samples}_{build}_.bed

Conda software stack

samtools ==1.10
bedtools =2.29.0

Code

1
2
3

        (samtools view {params.q_trim} -b {input} |
        bedtools bamtobed > {output}) 2>{log}

Jobs 1

Input files

results/{samples}/{samples}_cutadapt.fastq.gz
resources/ref_genomes/{build}/Bowtie2/genome_{build}.1.bt2

Output files

results/{samples}/{samples}_cutadapt_se_{build}.bam

Conda software stack

bowtie2 ==2.4.1
samtools ==1.10

Code

1 2	`(bowtie2 --threads {threads} {params.extra} -x {params.ref_genome} -U {input.sample[0]} \| samtools view -Sbh -o {output} -) 2>{log}`

Jobs 1

Input files

resources/samples/{samples}.fastq.gz

Output files

results/{samples}/{samples}_cutadapt.fastq.gz
results/{samples}/{samples}_cutadapt.qc.txt

Conda software stack

cutadapt ==2.10

Code

"""Snakemake wrapper for trimming single-end reads using cutadapt."""

__author__ = "Julian de Ruiter"
__copyright__ = "Copyright 2017, Julian de Ruiter"
__email__ = "julianderuiter@gmail.com"
__license__ = "MIT"


from snakemake.shell import shell


n = len(snakemake.input)
assert n == 1, "Input must contain 1 (single-end) element."

extra = snakemake.params.get("extra", "")
adapters = snakemake.params.get("adapters", "")
log = snakemake.log_fmt_shell(stdout=False, stderr=True)

assert (
    extra != "" or adapters != ""
), "No options provided to cutadapt. Please use 'params: adapters=' or 'params: extra='."

shell(
    "cutadapt"
    " {snakemake.params.adapters}"
    " {snakemake.params.extra}"
    " -j {snakemake.threads}"
    " -o {snakemake.output.fastq}"
    " {snakemake.input[0]}"
    " > {snakemake.output.qc} {log}"
)

Jobs 1

Input files

results/{samples}/{samples}_{build}_sorted_chr.bed

Output files

results/{samples}/{samples}_{build}_length_distribution.txt

Code

1
2
3

        awk '{{print $3-$2}}' {input} | sort -k1,1n | uniq -c | 
        sed 's/\s\s*/ /g' | awk '{{print $2"\t"$1}}' > {output}

Jobs 1

Input files

results/{samples}/{samples}_{build}_sim_dinucleotideTable.txt
results/{samples}/{samples}_{build}_sim_nucleotideTable.txt

Output files

results/{samples}/{samples}_{build}_sim_dinucleotideTable.png
results/{samples}/{samples}_{build}_sim_nucleotideTable.png

Conda software stack

conda-forge::r-base ==3.5.1
conda-forge::r-argparser ==0.4
conda-forge::r-ggplot2 ==3.0.0
conda-forge::r-dplyr ==0.8.0.1
conda-forge::r-tidyr ==0.8.2

Code

        Rscript workflow/scripts/nucleotidePlot.r         -i {input.dinuc}         -k 2         -s {params.name}         -f {params.motif}         -o {output.dinuc} 2>{log}

        Rscript workflow/scripts/nucleotidePlot.r         -i {input.nuc}         -k 1         -s {params.name}         -o {output.nuc} 2>>{log}

Jobs 1

Input files

results/{samples}/{samples}_{build}_xr_sim.fa

Output files

results/{samples}/{samples}_{build}_sim_dinucleotideTable.txt
results/{samples}/{samples}_{build}_sim_nucleotideTable.txt
results/{samples}/{samples}_{build}_xr_sim_filt.fa

Code

        len="$(awk 'BEGIN{{RS=">";ORS=""}}{{print length($2)"\n"}}' {input} | sort | uniq -c | sort -k1,1n | tail -1 | awk '{{print $2}}')"

        awk -v num="$len" 'BEGIN{{RS=">";ORS=""}}length($2)==num{{print ">"$0}}' {input} > {output.filt}

        workflow/scripts/fa2kmerAbundanceTable.py         -i {output.filt}         -k 2         -o {output.dinuc} 2>{log}

        workflow/scripts/fa2kmerAbundanceTable.py         -i {output.filt}         -k 1         -o {output.nuc} 2>>{log}

Jobs 1

Input files

results/{samples}/{samples}_{build}_sorted_chr.bed
resources/ref_genomes/{build}/genome_{build}.fa
resources/ref_genomes/{build}/genome_{build}.ron

Output files

results/{samples}/{samples}_{build}_sorted_chr.fa
results/{samples}/{samples}_{build}_xr_sim.fa

Conda software stack

bedtools =2.29.0

Code

        bedtools getfasta         -fi {input.genome}         -bed {input.bed}         -fo {output.fa}         -s

        boquila         --fasta {output.fa}         --ref {input.genome}         --regions {input.regions}         > {output.sim}

Jobs 1

Input files

resources/ref_genomes/{build}/genome_{build}.fa.fai

Output files

resources/ref_genomes/{build}/genome_{build}.ron

Code

1	`python3 workflow/scripts/idx2ron.py -i {input} -o {output}`

Jobs 1

Input files

results/{samples}/{samples}_{build}_length_distribution.txt

Output files

results/{samples}/{samples}_{build}_length_distribution.png

Conda software stack

conda-forge::r-base ==3.5.1
conda-forge::r-argparser ==0.4
conda-forge::r-ggplot2 ==3.0.0
conda-forge::r-dplyr ==0.8.0.1
conda-forge::r-tidyr ==0.8.2

Code

1 2	`Rscript workflow/scripts/lengthDistPlot.r -i {input} -s {params} -o {output} 2>{log}`

Jobs 2

Input files

results/{samples}/{samples}_{build}_sorted_xr_{strand}.bdg
resources/ref_genomes/{build}/genome_{build}.fa.fai

Output files

results/{samples}/{samples}_{build}_sorted_xr_{strand}.bw

Conda software stack

ucsc-bedgraphtobigwig == 377

Code

1 2	`bedGraphToBigWig {input.bdg} {input.index} {output}`

Jobs 1

Input files

results/{samples}/{samples}_{build}_sorted_plus.bed
results/{samples}/{samples}_{build}_sorted_minus.bed
resources/ref_genomes/{build}/genome_{build}.fa.fai

Output files

results/{samples}/{samples}_{build}_sorted_xr_plus.bdg
results/{samples}/{samples}_{build}_sorted_xr_minus.bdg

Conda software stack

bedtools =2.29.0

Code

        bedtools genomecov -i {input.plus} -g {input.ref_genome} -bg -scale         $(echo {params.read} | awk '{{print 1000000/$1}}') > {output.plus}

        bedtools genomecov -i {input.minus} -g {input.ref_genome} -bg -scale         $(echo {params.read} | awk '{{print 1000000/$1}}') > {output.minus}

Jobs 1

Input files

results/{samples}/{samples}_{build}_sorted_chr.bed

Output files

results/{samples}/{samples}_{build}_sorted_plus.bed
results/{samples}/{samples}_{build}_sorted_minus.bed

Code

1
2
3

        awk '{{if($6=="+"){{print}}}}' {input} > {output.plus}
        awk '{{if($6=="-"){{print}}}}' {input} > {output.minus}